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Abstract. The concept of molecular mechanics force held has been widely accepted nowa¬ 
days for studying various processes in biomolecular systems. In this paper, we suggest a 
modihcation for the standard CHARMM force held that permits simulations of systems 
with dynamically changing molecular topologies. The implementation of the modihed force 
held was carried out in the popular program MBN Explorer, and, to support the devel¬ 
opment, we provide several illustrative case studies where dynamical topology is necessary. 
In particular, it is shown that the modihed molecular mechanics force held can be ap¬ 
plied for studying processes where rupture of chemical bonds plays an essential role, e.g., 
in irradiation- or collision-induced damage, and also in transformation and fragmentation 


processes involving biomolecular systems. 

1 Introduction 

Nowadays, it has become feasible to study structure 
and dynamics of molecular systems that constitute 
of millions of atoms and evolve on time scales 
up to hundreds of nanoseconds Q by employing the 
classical molecular mechanics (MM) approach. In 
this approach, a molecular system is treated classi¬ 
cally, so that constituent atoms interact with each 
other through a parametric phenomenological po¬ 
tential that is governed by the type of individual 
atoms and by the network of chemical bonds be¬ 
tween them. This network defines a so-called molec¬ 
ular topology, that is a set of rules that impose 
constraints in the system and permit maintaining 
its natural shape, mechanical, and thermodynam¬ 
ical properties. The MM method has been widely 
used throughout the last decades and imple¬ 

mented in the well-established computational pack¬ 
ages, such as CHARMM 0, AMBER [|, GRO- 
MACS 0, and NAMD 0. 

MBN Explorer (www.mbnexplorer.com) m 
is an alternative, emerging software for the simu¬ 
lation of complex biomolecular, nano- and meso¬ 
scopic systems. It is suitable for classical molecular 
dynamics (MD), Monte Carlo [l2l-[lq and relativis¬ 
tic dynamics simulations |17H^ of a large range 


of molecular systems, such as nano- [U, and 
biological systems, nanostructured materials [2^, 
[^ . composite/hybrid materials Esl-I^. gases, liq¬ 
uids, solids and various interfaces [2^1^, with the 
sizes ranging from atomic to mesoscopic. Among 
other applications, MBN Explorer can be used 
to simulate thermo-mechanical damage of a bio¬ 
logical medium, e.g. a DNA nucleosome, which is 
caused by the propagation of a shock wave initi¬ 
ated by irradiation with fast ions M- The results 
of such simulations are used then to evaluate the 
efficiency of radiation with different projectiles (3^ 
within the framework of the multiscale approach 
to the physics of radiation damage and can be 
^applied in the field of ion-beam cancer therapy (3^ 

Despite numerous successes, the conventional 
MM method is primarily capable of studying pro¬ 
cesses where chemical reactions do not take place. 
This leads to significant limitations of the method 
and makes it practically unsuitable for studying 
highly non-equilibrium processes in biomolecular 
systems, e.g. thermo-mechanical biodamage. This 
particular example involves rupture and formation 
of covalent bonds that cannot be simulated by the 
conventional MM method due to a fixed topology 
of the system. 



Simulation of the rupture and formation of co¬ 
valent bonds can be performed by using Quantum 
Mechanical/Molecular Mechanical (QM/MM) me¬ 
thods or ab initio MD simulations (36l - l3q . Both 
methods are computationally rather demanding and, 
thus, the ab initio approach is used typically for 
studying fragmentation of small biomolecules, such 
as DNA nucleobases or nucleotides [H, The 
size of such systems is far from the typical sizes of 
systems of biological relevance, consisting of hun¬ 
dred thousands of atoms, and more. This problem 
is addressed to some extent in QM/MM methods 
where a core part of a large biomolecular system is 
described quantum mechanically while all the sur¬ 
roundings are described classically u sing, fo r exam¬ 
ple, the conventional MM method [dlTld^. Thus, 
the rupture or formation of covalent bonds can be 
simulated only in a small part of the system, which 
is treated quantum mechanically. 

In this paper, we present an extended version of 
the conventional MM method, which has been re¬ 
cently implemented in MBN Explorer, and de¬ 
monstrate that this extension describes correctly 
the dynamically changing molecular topology of a 
system within the classical MD framework. The 
presented modification takes into account additional 
parameters of the system, such as dissociation en¬ 
ergy of bonds, bonds multiplicity and the valence 
of atoms. The functional form of the interatomic 
interactions is also adjusted to account for the fi¬ 
nite dissociation energy of the chemical bonds. Fi¬ 
nally, three examples of simulations with the ex¬ 
tended MM method are presented to demonstrate 
a proof of principle for utilizing the force field with 
dynamic topology. The first two case studies illus¬ 
trate the processes of rupture and formation of co¬ 
valent bonds in a small biomolecule, namely an ala¬ 
nine dipeptide, which is one of the simplest building 
blocks of larger biomolecular systems like polypep¬ 
tides or proteins. Having proven the force field to 
work on a simple dipeptide allows us to generalize 
the framework towards macromolecules. This will 
allow for studying the systems of biologically rele¬ 
vant sizes, on the time scales which are not acces¬ 
sible by means of ab initio methods. The last ex¬ 
ample illustrates the process of water splitting and 
the evolution of chemical equilibrium. It is demon¬ 
strated that the results of the simulation are in a 
reasonable quantitative agreement with those of the 
analytical calculations. 


2 Theoretical approach 

2.1 Molecular mechanics potential 

The MM potential represents a phenomenological 
parametrization of the potential energy of a system 
and is widely used to describe structure and prop¬ 


erties of macromolecular systems, such as polypep - 
tides l43l- i^ pr oteins (ttI - I^ . DNA |32l! [5^453 |. 
lipids |51L I55l457j| . and many others. All physically 
important interactions in a system, i.e. both cova¬ 
lent and non-bonded long-range interactions, are 
accounted for using a simple parametric form, so 
that the total energy of the system reads as: 

Utot — Ucov T ffvdW T f^Coul ■ (1) 


The terms on the right-hand side describe the co¬ 
valent, van der Waals, and electrostatic Coulomb 
interactions, respectively. The van der Waals inter¬ 
action between two neutral atoms or molecules is 
usually modeled by the Lennard-Jones (LJ) poten¬ 
tial: 


C/Lj = e 



( 2 ) 


where e is the depth of the potential energy well, rg 
is the equilibrium distance, and is the distance 
between atoms. The function f/cov parameterizes 
the covalent interactions through a number of em¬ 
pirical parameters and fitting functions. MBN Ex¬ 
plorer implies the following parametrization for 
the Ucov term: 
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It is important to stress here that MM potential is 
only applicable to the systems with a predefined 
topology, i.e. a set of rules that define chemical 
bonds in the system. The energy of the covalent 
interactions is then modeled as a sum over all such 
interactions. The first and the second terms on the 
right-hand side of Eq. dJ) describe the components 
of the potential energy of the system arising due 
to stretching of the bonds between two atoms and 
due to variation of angles between every topologi¬ 
cally defined triplet of atoms, respectively (see Fig¬ 
ure [TJi). The third term is known as the torsion 
energy, which is characterized through a dihedral 
angle formed by every four atoms connected via co¬ 
valent chemical bonds. The last term describes the 
so-called improper dihedral angles that are used in 
the molecular topology to maintain planarity. 

A widely used MM potential, called CHARMM 
[S^ . represents a set of parameters for the sim¬ 
ulation of structure and dynamics of bio-/macro- 
molecular systems. This force field employs har¬ 
monic approximation for describing the interatomic 
interactions, thereby limiting its applicability to 
small deformations of the molecular system. The 
form of the potential energy functions, which are 


2 





Fig. 1. (a) Internal coordinates describing molecu¬ 
lar mechanics interactions: rij governs bond stretch¬ 
ing, Oijk represents the angular term, Xijki gives the 
dihedral angle, and the small out-of-plane angle Sijki is 
governed by the so-called improper dihedral angle, (b) 
Dependencies of the potential energy on coordinates 
used in the molecular mechanics potential, Eq. Q , de¬ 
scribing the bonded, angular, dihedral angular, and im¬ 
proper dihedral angular interactions. 


used to describe the components of Ucov in the stan¬ 
dard CHARMM force field, is demonstrated in Fig¬ 
ure [Br. In the case of substantial deformations, the 
interaction forces should decrease to zero as the va¬ 
lence bonds rupture. The rupture of valence bonds 
should also cause the involved angular and dihe¬ 
dral interactions to vanish as well. In this work, we 
report on advances in the MEN Explorer devel¬ 
opment, that permit classical MD simulations of 
the rupture of covalent bonds by using a dissocia¬ 
tive CHARMM potential. This approach goes be¬ 
yond the harmonic approximation, thus describing 
the physics of molecular dissociation more accu¬ 
rately, and permits construction of dynamic molec¬ 
ular topology, which instructs MEN Explorer 
how the existing covalent bonds can break and new 
covalent bonds can be formed. These features make 
MEN Explorer rather unique, for example, for 
simulating irradiation- and collision-induced bio¬ 
damage by means of classical MD. To the best of 
our knowledge, presently there is no similar ap¬ 


proach implemented in other software for classical 
MD simulations and, therefore, the implemented 
algorithms are seen as unique know-how of MEN 
Explorer. The adopted methodology is described 
in detail further in this section. 


2.2 Rupture of covalent bonds 

In order to model rupture of covalent bonds in the 
CHARMM force field, MEN Explorer uses mod¬ 
ified interaction potentials which describe the inter¬ 
actions of atoms connected by chemical bonds. The 
standard CHARMM force field treats the covalent 
interactions within the harmonic approximation as 

[/<;■“"» = (4) 

Here k]’, is the force constant of the bond stretch- 

‘‘J 

ing, rij is the distance between atoms i and j, and 
the parameter rg is the equilibrium covalent bond 
length. The above parametrization describes well 
the bond stretching regime in the case of small de¬ 
viations from rg but gives an erroneous result for 
larger distortions. For a satisfactory description of 
the covalent bond rupture it is reasonable to substi¬ 
tute the parabolic potential with the Morse poten¬ 
tial. It requires one additional parameter, as com¬ 
pared to the harmonic approximation Q , and takes 
into account the energy of bond dissociation. Po¬ 
tential energy of the system of two atoms interact¬ 
ing via the Morse potential reads as: 


Uuirij) = Di 


,-2Pij{rij-ro) _ 2g-/3i3(’’«-ro) 


( 5 ) 

where Dij is the bond dissociation energy and the 
parameter Pij determines steepness of the poten¬ 
tial. It follows from Eq. ([^ that C/M(rg) = —Dij- 


Let us consider a small deformation of the cova¬ 
lent bond from its equilibrium distance, — rg <C 
rg. In this case, the energy of the bond can be ap¬ 
proximated harmonically as: 


UM(rij) Ki -Dij -I- Pij'^D^j{rij - r^fk , (6) 


so that 

N = ■ ( 7 ) 

This expression defines Pij for a certain covalent 
bond and relates it to the value of k^j used in the 
standard CHARMM force field. 

Figure [2] illustrates the Morse potential which 
models the interaction of two carbon atoms, having 
the CN7 - CN8E type according to the CHARMM 
nomenclature. This covalent bond occurs, for exam¬ 
ple, between the C4’-C5’ atoms of ribose. At small 
deviations from rg, the Morse potential (red curve) 
and the harmonic approximation (green curve) are 
close to each other. With increasing interatomic 
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Fig. 2. The pairwise carbon-carbon (CHARMM type 
CN7 - CN8B) interaction potential in harmonic (1) and 
Morse approximations (2). The van der Waals interac¬ 
tion between the two atoms is illustrated with the blue 
line denoted as (3). 



Fig. 3. The switching function cr(rij) calculated for the 
carbon-carbon (CHARMM type CA -CA) interaction 
using Eq. uni. The function is used to rupture angular 
interactions in the modified CHARMM force field. 


distance (r^ rg), the atoms start to interact 
through polarization forces modeled through the 
Lennard-Jones potential. The comparison between 
the Morse and the Lennard-Jones potentials at larger 
distances is shown in the inset. It follows that both 
potentials are close to each other at the distances 
of about the van der Waals contact distance for the 
non-bonded interactions, which is about 2 A for the 
considered carbon atoms. 


2.3 Rupture of valence angles 

The rupture of chemical bonds in the course of sim¬ 
ulation automatically employs an improved poten¬ 
tial for the valence angles. In the CHARMM force 
field, the potential associated with the change of a 
valence angle between bonds with indices ij and j k 
reads as: 

= Kjk {0,jk - 0of , ( 8 ) 

where and 9g are parameters of the potential, 
and dijk is the actual value of the angle formed 
by the three atoms. This potential grows rapidly 
with increasing the angle, and it may lead to non¬ 
physical results when modeling the covalent bond 
rupture. In order to avoid such cases, in the mod¬ 
ified force field the harmonic potential (|H]) is sub¬ 
stituted with an alternative parametrization, which 
reads as: 

[1 - cos(%fe - 0o)] . (9) 

At small variations of the valence angle, this para¬ 
metrization is identical to the harmonic approxima¬ 
tion Q used in the standard CHARMM force field. 
For larger values of the angle, the new parametriza¬ 
tion dni) defines an energy threshold which becomes 
important for an accurate modeling of bond break¬ 
age. 


The rupture of a covalent bond is accompanied 
by a rupture of the angular interactions associated 
with this bond. The effect of bond breakage on the 
angular potential can be described through a spe¬ 
cial function a{rij), defined as 

^ [Pijirij - r*j)] } , ( 10 ) 

with r*j = + ro)/2. Figure [3] illustrates that 

has a form of a smoothed step function. This 
function introduces a correction to the angular in¬ 
teraction potential, assuming that the distance be¬ 
tween two atoms involved in an angular interaction 
increases from the equilibrium value, tq, up to the 
van der Waals contact value, . Since an an- 
gular interaction depends on two bonds connect¬ 
ing the atoms with indices ij and ik, the potential 
energy, describing the valence angular interaction 
that is subject to rupture, is parameterized as 

( 11 ) 

As seen from this expression, the angular poten¬ 
tial decreases with the increase of the bond length 
between any of the two pairs of atoms ij or jk. 
For exemplary purposes, in Figure 0^ we show the 
CN8B - ON2 - P angular potential, which arises, 
for instance, when modeling DNA nucleotides. The 
presented angular potential was calculated using 
Eq. (fTTll assuming the breakage of the bond be¬ 
tween the oxygen and the phosphorous atoms. For 
the sake of illustration, the CN8B - ON2 bond 
length in this case was taken equal to its equilib¬ 
rium value tq. 

2.4 Rupture of dihedral interactions 

“Dihedral” interactions arise in the conventional 
MM potential due to the change of the dihedral an¬ 
gles between every four topologically defined atoms. 
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Fig. 4. (a) The CN8B - ON2 - P angular potential 
calculated using Eq. m with account for the ON2 - P 
bond rupture, (b) The CN4 - P - ON2 - CN7 dihedral 
potential calculated using Eq. m with account for the 
ON2 - P bond rupture. 


Let us consider a quadruple of atoms with indices i, 
j, k and I (see Figure[T]), bound through an interac¬ 
tion which is governed by a change of the dihedral 
angle. In this case, the dihedral angle stands for 
the angle between the plane, formed by the atoms 
i, j and k, and the plane, formed by the atoms j, k 
and 1. In the harmonic approximation, the dihedral- 
energy contribution reads as: 

J^puhedral) ^ COs { n,,kl Xijkl - S , jkl )] , 

( 12 ) 

where Uijki and Sijki are parameters of the 

potential, and Xijki is the angle between the planes 
formed by atoms f, j, k and j, fc, 1. 

The dihedral interactions also become disturbed 
upon covalent bond rupture; therefore, Eq. m 
should be modified to properly account for this ef¬ 
fect. The rupture of a dihedral interaction between 
a quadruple of atoms i, j, k and I should take into 
account three bonds that contribute to this inter¬ 
action. Thus, the potential energy describing the 
dihedral interaction with account for the bond rup¬ 
ture reads as: 

cS"”" = »(r«) »(d0 »(r«) , (13) 

where is the potential (IT^ describing 

the dihedral interaction within the framework of 


the standard CHARMM force field. The functions 
cr{rij), cr{rjk), and a{rki) are defined by Eq. (dH]); 
they are used to limit the dihedral interaction upon 
increasing the corresponding bond length. FigureHlD 
shows a typical profile of a dihedral potential with 
accounting for the bond rupture. In this case, we 
have considered the CN4 - P - ON2 - CN7 di¬ 
hedral interaction where the middle ON2 - P bond 
was broken. This interaction is also important when 
modeling bond breakages in DNA nucleotides. 


2.5 Formation of new bonds 

Rupture of covalent bonds leads to formation of in¬ 
dividual atoms, radicals or smaller molecular frag¬ 
ments. In order to properly simulate the chemical 
balance in the system, one should allow for the for¬ 
mation of new bonds. 

In MEN Explorer, after the rupture of a chem¬ 
ical bond between two atoms, these atoms are placed 
in a special list of chemically active atoms. Only the 
atoms from this list can participate in chemical re¬ 
actions and form new bonds. For each atom from 
this list, the number of possible molecular bonds is 
stored and determined by its valence. 

In order to create new bonds in the system, the 
list of chemically active atoms is evaluated at each 
simulation step, and the neighboring atoms are se¬ 
lected. A chemical bond is formed between a pair of 
atoms provided that the following conditions have 
been met: (i) the parameters of the bond, e.g., the 
equilibrium distance and the value of bond forma¬ 
tion energy for this combination of atoms are de¬ 
fined in the simulation input, (ii) atoms are mod¬ 
eled as bound through the Morse potential, and 
(iii) the distance between the atoms is less than 
the predefined cutoff/capture radius, which is an 
independent parameter for each bond type used to 
speed-up the simulations. If all these conditions are 
met simultaneously in the simulation, the bond is 
created and the system’s topology is updated. For 
each new covalent bond, the neighboring atoms are 
analyzed. If the parameters of angular bond are 
defined for some group of atoms such bond should 
also be formed. 


2.6 Redistribution of partial charges 

The rupture and formation of bonds leads to re¬ 
configuration of molecular systems and to redistri¬ 
bution of partial charges of atoms. The modifica¬ 
tion of MEN Explorer presented in this paper 
accounts for this phenomenon. The charge redis¬ 
tribution should obey the following conditions: the 
total charge of the systems is conserved, the to¬ 
tal charge of each individual molecule is an integer 
in atomic system of units, i.e. an integer number of 
the electron charge. MEN Explorer supports the 
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Fig. 5. Snapshots illustrating dynamics of alanine dipeptide and the C-N bond rupture simulated with the 
harmonic (left) and Morse (right) potentials at 0 fs (top), 8 fs (middle) and 20 fs (bottom). 


two options for modeling the charge redistribution 
process: (i) a general (default) one applicable to 
any molecule, and (ii) a special one where charges 
within molecules are redistributed according to the 
known electronic configurations. 


The default mechanism of charge redistribution 
is activated upon rupture or formation of cova¬ 
lent bonds. In the case of a bond rupture, two 
newly emerged fragments of a molecular system 
have likely non-compensated, non-integer charges. 
The total charge of each of the newly emerged frag¬ 
ment is thus rounded to the closest integer value 
and the charge difference is transferred from one 
fragment to another. This difference is redistributed 
evenly among all atoms of the fragments. Upon the 
formation of new a bond, the charge is redistributed 
inside the newly created molecule in order to lower 
the values of partial charges preserving the initial 
sum of charges. 


3 Numerical Results 

Let us now consider the case studies that go beyond 
the standard MM methodology. The first example 
illustrates the rupture of a single C-N bond in an 
alanine dipeptide molecule. In the second example, 
we simulate the reverse process of the new bond 
formation. These two case studies demonstrate a 
proof of principle for the modified MM force field to 
describe the dynamically changing molecular topol¬ 
ogy of the system within the classical MD frame¬ 
work. The third example is devoted to the inves¬ 
tigation of the process of water splitting at high 
temperatures. In this case, the breakage of chem¬ 
ical bonds and the formation of new ones lead to 
establishing the chemical equilibrium in the system. 


3.1 Fragmentation of alanine dipeptide 

To illustrate the bond breakage, we have simulated 
the dynamics of alanine dipeptide consisting of 20 
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atoms, solvated in a simulation box with 95 wa¬ 
ter molecules. The alanine dipeptide molecule was 
considered with neutral terminals. 

In order to clearly illustrate the difference be¬ 
tween the standard CHARMM force field, utilizing 
the harmonic interatomic potential, and the disso¬ 
ciative CHARMM potential implemented in MEN 
Explorer, two simulations were carried out. In 
these simulations, the rupture of the central C-N 
bond in the dipeptide, leading to the formation of 
two isolated alanine molecules, was monitored. To 
facilitate the process, we have set the initial ve¬ 
locity of the C and N atoms high, which corre¬ 
sponds to an energy fluctuation sufficient for the 
bond rupture. In the simulation performed with 
the standard force field, the peptide bond is mod¬ 
eled through the harmonic potential, therefore, the 
bond cannot break. The behavior of the C-N bond 
in the harmonic approximation is illustrated in the 
left part of Figure [SJ and the corresponding atoms 
are marked with red circles. In this case, the dis¬ 
tance between the atoms oscillates around the equi¬ 
librium value as the atoms always return to their 
equilibrium positions. 



Fig. 6. Dependence of the C-N interactomic distance 
in alanine dipeptide as a function of the simulation 
time. 


In the second simulation, the Morse potential ([S]) 
was used for the description of the peptide bond. 
In this case the C and N atoms do not oscillate 
around an equilibrium position, and the structure 
of the system after 20 fs of simulation changes sig¬ 
nificantly from the one considered above (see the 
right part of Figure [S]). It is evident from the snap¬ 
shots that the distance between the atoms increases 
already after 8 fs. 

When the distance between the atoms exceeds 
a given cutoff radius (which is equal to 2.5 A in 
this example), the bond is considered as broken. 
Once this has happened, the carbon and the nitro¬ 
gen atoms remain interacting only via the electro¬ 
static potential and the van der Waals interactions. 


so that the two alanine molecules can diffuse apart. 
The charge redistribution does not happen in this 
case because both new fragments of the dipeptide 
were initially neutral. 

FigureElshows the interatomic distance between 
the carbon and the nitrogen atoms as a function of 
the simulation time. The equilibrium distance be¬ 
tween the atoms is = 1.354 A (dashed line). 

The figure demonstrates that in the case of the sim¬ 
ulation with the Morse potential, the interatomic 
distance monotonically increases indicating that the 
bond is broken and that two isolated alanine molecu¬ 
les drift apart. 


3.2 Binding of two alanine amino acids 





Fig. 7. (a) Two alanine molecules approaching each 
other to form a new C-N bond, (b) Dependence of the 
distance between C and N atoms for the two alanine 
molecules. 


The second example illustrates the process of 
binding two alanine molecules together into a single 
dipeptide through the formation of a new covalent 
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bond in the molecular system. In this case study, 
six isolated alanine amino acids surrounded by 54 
water molecules were placed in a small simulation 
box of 24 X 24 X 24 with periodic boundary con¬ 
ditions, and the dynamics of the system was simu¬ 
lated for 80 ps at a fixed temperature of 1000 K con¬ 
trolled by the Langevin thermostat with the damp¬ 
ing time constant of 1 fs. Each alanine molecule was 
modeled with unsaturated N- and C- termini, i.e. 
having two unpaired chemical bonds. 

In the course of the simulation, all distances be¬ 
tween the different termini of alanines were moni¬ 
tored. When the distance between a pair of termi¬ 
nal atoms became smaller than the predefined cut¬ 
off radius (equal to 3 A in this example), a new co¬ 
valent bond was considered to be formed. Figure[7^ 
illustrates a spacial conformation of two amino acids 
in the simulation leading to the formation of a new 
bond. Figure [7 }d gives the dependence of the dis¬ 
tance between C and N atoms for the two molecules 
shown in the upper part. At some point, this dis¬ 
tance becomes smaller than the cutoff radius (blue 
dashed line), and the two molecules become con¬ 
nected. Note that after 40 ps the distance between 
the C and N atoms oscillates around a constant 
value corresponding to the C-N bond equilibrium 
length. Since six alanines are considered in this 
simulation, more of them could self-assemble in a 
polypeptide chain but this would require longer 
simulation. In this system, each initial alanine mole¬ 
cule has a total charge equal to zero. Therefore, 
after the formation of a new molecule the charge 
redistribution step was not necessary. 


3.3 Water splitting 

The third case study concerns the process of wa¬ 
ter splitting. At elevated temperatures, water mol¬ 
ecules can dissociate forming different molecular 
products. In such a system, the following reactions 
take place: 


H 2 O -h 

(14) 

OH ^ 0"® -f H+® 

(15) 

H 2 H -H H 

(16) 

O 2 O-kO , 

(17) 


which describe both the mechanisms of rupture and 
formation of new chemical bonds. The superscripts 
+qi and —qi (i = 1,2) indicate that the reaction 
products in Eqs. (O and m carry some non-zero 
charge. The principle of detailed balance states that 
for a reversible process at the equilibrium a direct 
process should be equilibrated by its reverse pro¬ 
cess. This forms the equilibrium concentrations of 
reacting species which depend on the temperature 
of the medium. In order to determine the equilib¬ 
rium concentrations of molecules in our system. 


one should use the equations for the equilibrium 
concentrations of different molecular species [h^ . 
which are based on the Dalton’s law and the law of 
mass action [^, : 


Ah^oh.P 

AiHaO 

XkXqP 

^OH 

X 112 


= exp 


= exp 


= exp 


3L 

X 02 


exp 



(18) 

(19) 

( 20 ) 
( 21 ) 


where Ah, Aqh, Xq and AhjO are the relative 
concentrations of atoms and molecules of different 
types noted by the subscripts, AGi are Gibbs en¬ 
ergy values - the dissociation energies of the chem¬ 
ical reactions in Eqs. (Hi-dni), P is the total pres¬ 
sure in the system, R is the universal gas constant, 
and T is the temperature. Taking into account that 
the total number of hydrogen atoms in the studied 
system is twice the number of oxygen atoms, one 
obtains: 


2Ah2 + Ah -I- Aqh + 2Ah20 
= 2(2Ao2+Aoh + Ao + Ah2o) • (22) 

The normalization condition for the partial concen¬ 
trations/densities can be written accordingly as: 

Ah 2 + A02 + Ah 

+Aoh + Ao -I- AH 2 O = 1 ■ (23) 

In this case study, we have simulated water split¬ 
ting in a box of 130 x 130 x 130 A^ with peri¬ 
odic boundary conditions. The studied system con¬ 
sists of 10917 water molecules randomly placed in¬ 
side this box. Simulations for different tempera¬ 
tures were carried out using the Langevin thermo¬ 
stat with the damping time of 10 fs. The total sim¬ 
ulation time was set to 500 ps. We note that ac¬ 
counting for the dissociation and formation of new 
bonds in the MM method with dynamic topology 
does not affect significantly the performance of the 
simulation. The computational cost of the calcu¬ 
lation is about the same as for the conventional 
MD simulation with the standard CHARMM force 
field; the simulation of water splitting was carried 
out using a 12-core AMD workstation in 3 days. 
Due to the constant volume and density, and differ¬ 
ent temperatures at each simulation, the pressure 
in the box was changing from 2500 to 4000 bar 
which corresponds to the water vapor state with a 
high density. The simulated values of pressure and 
temperature were used in the analytic expressions, 
Eqs. (HU)-(1231), for estimating the equilibrium con¬ 
centrations of the reaction products. The values of 
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Fig. 8. Time evolution of the relative concentration of atoms Xn, Xq and molecules Xu^, Xoh, ^H 2 0 in the 
simulated system at T = 4800 K. Results of the simulation are shown with solid lines, dots show corresponding 
analytic values of concentrations. 
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Fig. 9. Comparison of the simulated dependence of 
the relative concentration of water molecules, XhjO, 
on temperature with the one obtained from the equi¬ 
librium analysis, Eqs. ([T5l)-(Pl). 




• 1^H20 simulation 
- -^H20 equilibrium 




Gi = G 2 = 117.9 kcal/mol, G 3 = 104 kcal/mol, 
G 4 = 119 kcal/mol were used. 

For each reaction (Hl-dlSl), the distribution of 
partial charges qi was calculated separately and set 
explicitly for all products. For water, the partial 
charge of oxygen and hydrogen atoms was set to 
—0.834 and 0.417, respectively, for a OH molecule 
the values of —0.375 and 0.375 were assumed. H 2 
and O 2 molecules in reactions (m and (HZl) were 
considered as neutral with zero partial charges. 

Figure[ 8 ]shows the time evolution of the relative 
concentration of atoms and molecules in the system 
at T = 4800 K. Solid lines present the results of 
the simulations, while symbols show the results of 
numerical solution of Eqs. m-m- 


Initially, all molecules in the system are water 
molecules. Frequent random collisions lead to their 
dissociation into H and OH. During the first 10 ps, 
these are the two main reaction products and their 
relative concentration coincide. After 10 ps the con¬ 
centration of hydrogen atoms becomes significant 
and the reaction H -|- H H 2 becomes probable. 
At this instance the concentration of H 2 increases 
significantly and the concentration of H drops. Af¬ 
ter 100 ps of simulation, the concentrations of H, 
H 2 , OH and H 2 O become nearly constant and close 
to the equilibrium values obtained from Eqs. m- 
(E31). However due to the low concentration of oxy¬ 
gen atoms the formation of O 2 molecules is a rel¬ 
atively slow process and requires significant time 
which we could not reach in our simulation. This 
process starts after about 100 ps and does not reach 
the saturation stage by the end of the simulation. 

To further characterize the interconversion of 
chemical reactions in the system, in Eigure [5] we 
present the comparison of the simulated temper¬ 
ature dependence of the relative concentration of 
water molecules, A'h 20 ; with the one obtained from 
the equilibrium analysis, Eqs. (1151) - (1251) . The in¬ 
crease of temperature leads to the splitting of water 
molecules and their number decreases. The com¬ 
parison shows that the two curves show practically 
identical behavior. Small deviations can be attribu¬ 
ted to the limited simulation time and to the high 
pressure in the system, as both of these factors lead 
to the decrease of the number of split molecules. 

A similar analysis was performed for Ah , Ah 2 , 
Aoh, and Xq. Figure [15] illustrates the dependence 
of relative concentration on temperature for the dif¬ 
ferent products of dissociation. These dependencies 
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Fig. 10. Comparison of simulated dependencies of Xh, 
Xh 2 , Xoh, Xo on temperature with those obtained 
from the equilibrium analysis, Eqs. (Il 8 j-(l 23 ll. Dashed 
lines correspond to the results of the equilibrium anal¬ 
ysis, solid curves correspond to the results of simula¬ 
tions. 


are compared with those obtained from the equilib¬ 
rium analysis based on Eqs. dm)-®. The figure 
demonstrates a reasonable agreement of the two ap¬ 
proaches, especially at high temperatures where the 
simulated values Xu nearly coincide with those ob¬ 
tained from the equilibrium analysis, and the values 
of Xh 2 and Xqh are rather close. For smaller tem¬ 
peratures, the deviation of the curves increases, pri¬ 
marily because of the insufficient simulation time. 
At higher temperatures, the molecules diffuse faster, 
thus allowing to achieve sooner the chemical equi¬ 
librium in the system. 


4 Conclusions 

This paper reports on an important development 
of the MEN Explorer software package, namely 
on the implementation of reactive force fields based 
on the modified CHARMM force field. This mod¬ 
ification allows one to simulate and analyze chem¬ 
ical reactions and transformations involving cova¬ 
lent bond breakages in bio- and organic molecules 
and molecular systems [^, |^. This is obviously 
a very large research domain and we believe that 
the implementations reported in this work will open 
many new possibilities for the research analysis bas¬ 
ed on the computational modeling of the mentioned 
systems. In this paper, we have not tried to explore 
all the possible applications of the new computa¬ 
tional tool, but instead focused on the most char¬ 
acteristic and illustrative examples. 

This work can be extended further in many dif¬ 
ferent ways. The introduced reactive force fields can 
be applied systematically to different systems and 
the results of simulations can be compared with the 
appropriate experiments [b^. The simulation 


results can also be validated through the compari¬ 
son with the corresponding results obtained through 
the quantum MD simulations of relatively small 
molecular systems [S^, . These and many more 

other possible studies based on the reported key im¬ 
plementation are the subject for further research. 
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